##########################################################
### Disorganized Political Violence:                   ###
### A Demonstration Case of Temperature and Insurgency ###
### Replication Code:                                  ###
### Temperature Response Function Results:             ###
### Mean Temperatures Robustness Check                 ###
### Andrew Shaver, Alex Bollfrass                      ###
### Date: 07/04/2022                                   ###
##########################################################

# Note: these results were generated using:
# 1) MacBook Pro (16-inch, 2021); Chip: Apple M1 Max; Memory: 64 GB; running
# macOS Monterey.
# 2) R 4.1.2 GUI 1.77 High Sierra build (8007)
# 3) RStudio 2022.07.1+554 "Spotted Wakerobin" Release (7872775ebddc40635780ca1ed238934c3345c5de, 2022-07-22) for macOS
# Mozilla/5.0 (Macintosh; Intel Mac OS X 12_1_0) AppleWebKit/537.36 (KHTML, like Gecko) QtWebEngine/5.12.10 Chrome/69.0.3497.128 Safari/537.36
# (Please see below for specific package versions.)

# Load required packages:
library(arm)
library(broom)
library(dplyr)
library(ggpubr)
library(ggplot2)
library(lfe)
library(lubridate)
library(MASS)
library(mgcv)
library(stargazer)
library(sandwich)
library(nnet)
library(scales)
library(gridExtra)
library(grid)

packageVersion("arm")
# ‘1.12.2’
packageVersion("broom")
# ‘0.7.12’
packageVersion("dplyr")
# ‘1.0.9’
packageVersion("ggpubr")
# ‘0.4.0’
packageVersion("ggplot2")
# ‘3.3.6’
packageVersion("lfe")
# ‘2.8.7.1’
packageVersion("lubridate")
# ‘1.8.0’
packageVersion("MASS")
# ‘7.3.55’
packageVersion("mgcv")
# ‘1.8.39’
packageVersion("stargazer")
# ‘5.2.3’
packageVersion("sandwich")
# ‘3.0.1’
packageVersion("nnet")
# ‘7.3.17’
packageVersion("scales")
# ‘1.2.0’
packageVersion("gridExtra")
# ‘2.3’
packageVersion("grid")
# ‘4.1.2’

# Load data:
# For primary temperature-violence analysis:
# Daily analysis:
iq_sig1 <- readRDS("~/Library/CloudStorage/Box-Box/TurningUpTheHeat1/IO_Replication_Materials/Processed Data/iraq_baghdad_basrah_processed_data.rds")
af_sig1 <- readRDS("~/Library/CloudStorage/Box-Box/TurningUpTheHeat1/IO_Replication_Materials/Processed Data/afghanstan_viol_districts_processed_data.rds")
# For attitude analysis:
bg_sur <- readRDS("~/Library/CloudStorage/Box-Box/TurningUpTheHeat1/IO_Replication_Materials/Processed Data/baghdad_survey.rds")
# For median analysis:
iq_median_temp <- readRDS("~/Library/CloudStorage/Box-Box/TurningUpTheHeat1/Affective_Motivations/MajorRevision/R3_Submission/ProcessedData/iq_daily_median.rds")

# Construct range of all temperatures observed across Afghanistan and Iraq:
all_observed_temps <- c(iq_sig1$temp, af_sig1$temp)
all_observed_temps <- as.data.frame(all_observed_temps)
colnames(all_observed_temps) <- "temperature"

# For temperature deviations for the Baghdad survey data, given that the unit 
# is not the day, match deviations from panel:
# Columns to keep:
to_keep <- c("date", "temp_1", "temp_2", "temp_3", "temp_4", 
             "temp_5", "temp_6", "temp_7", "daylight")
bg_sur <- plyr::join(bg_sur, iq_sig1[iq_sig1$governorate %in% "Baghdad", to_keep], by = "date")

# Create survey dataset removing top income earners:
bg_sur.s3.l <- bg_sur[bg_sur$income.w < 75000, ]

# Match median temperature data to the Iraq panel:
colnames(iq_median_temp)[2] <- "governorate"
iq_sig1 <- plyr::join(iq_sig1, iq_median_temp, by = c("date", "governorate"))

# Find minimum and maximum temperatures:
temperature_range <- round(min(iq_sig1$temp, af_sig1$temp, bg_sur$temp)):round(max(iq_sig1$temp, af_sig1$temp, bg_sur$temp))

###
### IRAQ.
###

# Baghdad and Basrah:
# Unconstrained Violence:
# OLS model results:

IQ.w.u.l <- felm(log(unconstrained+1) ~ temp + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + ied + constrained +
                   temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                   constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                   unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                   idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                   ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                   week + governorate | 0 | 0, data = iq_sig1)
IQ.w.u.q <- felm(log(unconstrained+1) ~ temp + temp2 + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + ied + constrained +
                   temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                   constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                   unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                   idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                   ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                   week + governorate | 0 | 0, data = iq_sig1)
IQ.w.u.q.lm <- lm(log(unconstrained+1) ~ temp + temp2 + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + ied + constrained +
                    temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                    constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                    unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                    idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                    ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                    as.factor(week) + as.factor(governorate), data = iq_sig1)
IQ.w.u.c <- felm(log(unconstrained+1) ~ temp + temp2 + temp3 + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + ied + constrained +
                   temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                   constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                   unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                   idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                   ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                   week + governorate | 0 | 0, data = iq_sig1)

IQ.m.u.l <- felm(log(unconstrained+1) ~ temp + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + ied + constrained +
                   temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                   constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                   unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                   idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                   ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                   month + governorate | 0 | 0, data = iq_sig1)
IQ.m.u.q <- felm(log(unconstrained+1) ~ temp + temp2 + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + ied + constrained +
                   temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                   constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                   unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                   idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                   ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                   month + governorate | 0 | 0, data = iq_sig1)
IQ.m.u.c <- felm(log(unconstrained+1) ~ temp + temp2 + temp3 + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + ied + constrained +
                   temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                   constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                   unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                   idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                   ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                   month + governorate | 0 | 0, data = iq_sig1)

# Constrained Violence:
IQ.w.c.l <- felm(log(constrained+1) ~ temp + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + ied + unconstrained +
                   temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                   constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                   unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                   idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                   ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                   week + governorate | 0 | 0, data = iq_sig1)
IQ.w.c.q <- felm(log(constrained+1) ~ temp + temp2 + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + ied + unconstrained +
                   temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                   constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                   unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                   idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                   ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                   week + governorate | 0 | 0, data = iq_sig1)
IQ.w.c.c <- felm(log(constrained+1) ~ temp + temp2 + temp3 + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + ied + unconstrained +
                   temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                   constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                   unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                   idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                   ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                   week + governorate | 0 | 0, data = iq_sig1)

IQ.m.c.l <- felm(log(constrained+1) ~ temp + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + ied + unconstrained +
                   temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                   constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                   unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                   idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                   ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                   month + governorate | 0 | 0, data = iq_sig1)
IQ.m.c.q <- felm(log(constrained+1) ~ temp + temp2 + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + ied + unconstrained +
                   temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                   constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                   unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                   idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                   ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                   month + governorate | 0 | 0, data = iq_sig1)
IQ.m.c.c <- felm(log(constrained+1) ~ temp + temp2 + temp3 + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + ied + unconstrained +
                   temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                   constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                   unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                   idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                   ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                   month + governorate | 0 | 0, data = iq_sig1)

# IED:
IQ.w.i.l <- felm(log(ied+1) ~ temp + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + constrained + unconstrained +
                   temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                   constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                   unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                   idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                   ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                   week + governorate | 0 | 0, data = iq_sig1)
IQ.w.i.q <- felm(log(ied+1) ~ temp + temp2 + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + constrained + unconstrained +
                   temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                   constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                   unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                   idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                   ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                   week + governorate | 0 | 0, data = iq_sig1)
IQ.w.i.c <- felm(log(ied+1) ~ temp + temp2 + temp3 + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + constrained + unconstrained +
                   temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                   constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                   unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                   idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                   ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                   week + governorate | 0 | 0, data = iq_sig1)

IQ.m.i.l <- felm(log(ied+1) ~ temp + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + constrained + unconstrained +
                   temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                   constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                   unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                   idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                   ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                   month + governorate | 0 | 0, data = iq_sig1)
IQ.m.i.q <- felm(log(ied+1) ~ temp + temp2 + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + constrained + unconstrained +
                   temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                   constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                   unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                   idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                   ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                   month + governorate | 0 | 0, data = iq_sig1)
IQ.m.i.c <- felm(log(ied+1) ~ temp + temp2 + temp3 + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + constrained + unconstrained +
                   temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                   constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                   unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                   idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                   ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                   month + governorate | 0 | 0, data = iq_sig1)

###
### AFGHANISTAN.
###

AF.w.u.l <- felm(log(df+1) ~ temp + prcp + wdsp + dewp + visib + mxspd + daylight + idf + ied +
                   temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                   df_1 + df_2 + df_3 + df_4 + df_5 + df_6 + df_7 +
                   idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                   ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                   week + station | 0 | month, data = af_sig1)
AF.w.u.q <- felm(log(df+1) ~ temp + temp2 + prcp + wdsp + dewp + visib + mxspd + daylight + idf + ied +
                   temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                   df_1 + df_2 + df_3 + df_4 + df_5 + df_6 + df_7 +
                   idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                   ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                   week + station | 0 | month, data = af_sig1)
AF.w.u.c <- felm(log(df+1) ~ temp + temp2 + temp3 + prcp + wdsp + dewp + visib + mxspd + daylight + idf + ied +
                   temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                   df_1 + df_2 + df_3 + df_4 + df_5 + df_6 + df_7 +
                   idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                   ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                   week + station | 0 | month, data = af_sig1)

AF.m.u.l <- felm(log(df+1) ~ temp + prcp + wdsp + dewp + visib + mxspd + daylight + idf + ied +
                   temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                   df_1 + df_2 + df_3 + df_4 + df_5 + df_6 + df_7 +
                   idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                   ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                   month + station | 0 | month, data = af_sig1)
AF.m.u.q <- felm(log(df+1) ~ temp + temp2 + prcp + wdsp + dewp + visib + mxspd + daylight + idf + ied +
                   temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                   df_1 + df_2 + df_3 + df_4 + df_5 + df_6 + df_7 +
                   idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                   ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                   month + station | 0 | month, data = af_sig1)
AF.m.u.c <- felm(log(df+1) ~ temp + temp2 + temp3 + prcp + wdsp + dewp + visib + mxspd + daylight + idf + ied +
                   temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                   df_1 + df_2 + df_3 + df_4 + df_5 + df_6 + df_7 +
                   idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                   ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                   month + station | 0 | month, data = af_sig1)

AF.w.i.l <- felm(log(ied+1) ~ temp + prcp + wdsp + dewp + visib + mxspd + daylight + df + idf +
                   temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                   df_1 + df_2 + df_3 + df_4 + df_5 + df_6 + df_7 +
                   idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                   ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                   week + station | 0 | month, data = af_sig1)
AF.w.i.q <- felm(log(ied+1) ~ temp + temp2 + prcp + wdsp + dewp + visib + mxspd + daylight + df + idf +
                   temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                   df_1 + df_2 + df_3 + df_4 + df_5 + df_6 + df_7 +
                   idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                   ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                   week + station | 0 | month, data = af_sig1)
AF.w.i.c <- felm(log(ied+1) ~ temp + temp2 + temp3 + prcp + wdsp + dewp + visib + mxspd + daylight + df + idf +
                   temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                   df_1 + df_2 + df_3 + df_4 + df_5 + df_6 + df_7 +
                   idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                   ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                   week + station | 0 | month, data = af_sig1)

AF.m.i.l <- felm(log(ied+1) ~ temp + prcp + wdsp + dewp + visib + mxspd + daylight + df + idf +
                   temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                   df_1 + df_2 + df_3 + df_4 + df_5 + df_6 + df_7 +
                   idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                   ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                   month + station | 0 | month, data = af_sig1)
AF.m.i.q <- felm(log(ied+1) ~ temp + temp2 + prcp + wdsp + dewp + visib + mxspd + daylight + df + idf +
                   temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                   df_1 + df_2 + df_3 + df_4 + df_5 + df_6 + df_7 +
                   idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                   ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                   month + station | 0 | month, data = af_sig1)
AF.m.i.c <- felm(log(ied+1) ~ temp + temp2 + temp3 + prcp + wdsp + dewp + visib + mxspd + daylight + df + idf +
                   temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                   df_1 + df_2 + df_3 + df_4 + df_5 + df_6 + df_7 +
                   idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                   ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                   month + station | 0 | month, data = af_sig1)

###
### Create dataset for prediction:
###

new_data_iq <- data.frame(temp = temperature_range, 
                          temp2 = (temperature_range)^2,
                          temp3 = (temperature_range)^3,
                          
                          prcp = mean(iq_sig1$prcp, na.rm = T), 
                          wdsp = mean(iq_sig1$wdsp, na.rm = T), 
                          dewp = mean(iq_sig1$dewp, na.rm = T), 
                          visib = mean(iq_sig1$visib, na.rm = T), 
                          mxspd = mean(iq_sig1$mxspd, na.rm = T), 
                          daylight = mean(iq_sig1$daylight, na.rm = T), 
                          hoursofpower = mean(iq_sig1$hoursofpower, na.rm = T), 
                          
                          unconstrained = mean(iq_sig1$unconstrained), 
                          constrained = mean(iq_sig1$constrained),  
                          idf = mean(iq_sig1$idf),  
                          ied = mean(iq_sig1$ied),  
                          
                          temp_1 = mean(iq_sig1$temp_1, na.rm = T),
                          temp_2 = mean(iq_sig1$temp_2, na.rm = T),
                          temp_3 = mean(iq_sig1$temp_3, na.rm = T),
                          temp_4 = mean(iq_sig1$temp_4, na.rm = T),
                          temp_5 = mean(iq_sig1$temp_5, na.rm = T),
                          temp_6 = mean(iq_sig1$temp_6, na.rm = T),
                          temp_7 = mean(iq_sig1$temp_7, na.rm = T),
                          
                          unconstrained_1 = mean(iq_sig1$unconstrained_1, na.rm = T),
                          unconstrained_2 = mean(iq_sig1$unconstrained_2, na.rm = T),
                          unconstrained_3 = mean(iq_sig1$unconstrained_3, na.rm = T),
                          unconstrained_4 = mean(iq_sig1$unconstrained_4, na.rm = T),
                          unconstrained_5 = mean(iq_sig1$unconstrained_5, na.rm = T),
                          unconstrained_6 = mean(iq_sig1$unconstrained_6, na.rm = T),
                          unconstrained_7 = mean(iq_sig1$unconstrained_7, na.rm = T),
                          
                          constrained_1 = mean(iq_sig1$constrained_1, na.rm = T),
                          constrained_2 = mean(iq_sig1$constrained_2, na.rm = T),
                          constrained_3 = mean(iq_sig1$constrained_3, na.rm = T),
                          constrained_4 = mean(iq_sig1$constrained_4, na.rm = T),
                          constrained_5 = mean(iq_sig1$constrained_5, na.rm = T),
                          constrained_6 = mean(iq_sig1$constrained_6, na.rm = T),
                          constrained_7 = mean(iq_sig1$constrained_7, na.rm = T),
                          
                          idf_1 = mean(iq_sig1$idf_1, na.rm = T),
                          idf_2 = mean(iq_sig1$idf_2, na.rm = T),
                          idf_3 = mean(iq_sig1$idf_3, na.rm = T),
                          idf_4 = mean(iq_sig1$idf_4, na.rm = T),
                          idf_5 = mean(iq_sig1$idf_5, na.rm = T),
                          idf_6 = mean(iq_sig1$idf_6, na.rm = T),
                          idf_7 = mean(iq_sig1$idf_7, na.rm = T),
                          
                          ied_1 = mean(iq_sig1$ied_1, na.rm = T),
                          ied_2 = mean(iq_sig1$ied_2, na.rm = T),
                          ied_3 = mean(iq_sig1$ied_3, na.rm = T),
                          ied_4 = mean(iq_sig1$ied_4, na.rm = T),
                          ied_5 = mean(iq_sig1$ied_5, na.rm = T),
                          ied_6 = mean(iq_sig1$ied_6, na.rm = T),
                          ied_7 = mean(iq_sig1$ied_7, na.rm = T),
                          
                          week = median(iq_sig1$week),
                          month = median(iq_sig1$month),
                          governorate = "Baghdad")

new_data_af <- data.frame(temp = temperature_range, 
                          temp2 = (temperature_range)^2,
                          temp3 = (temperature_range)^3,
                          
                          prcp = mean(af_sig1$prcp, na.rm = T), 
                          wdsp = mean(af_sig1$wdsp, na.rm = T), 
                          dewp = mean(af_sig1$dewp, na.rm = T), 
                          visib = mean(af_sig1$visib, na.rm = T), 
                          mxspd = mean(af_sig1$mxspd, na.rm = T), 
                          daylight = mean(af_sig1$daylight, na.rm = T), 
                          
                          df = mean(af_sig1$df), 
                          idf = mean(af_sig1$idf),  
                          ied = mean(af_sig1$ied), 
                          
                          temp_1 = mean(af_sig1$temp_1, na.rm = T),
                          temp_2 = mean(af_sig1$temp_2, na.rm = T),
                          temp_3 = mean(af_sig1$temp_3, na.rm = T),
                          temp_4 = mean(af_sig1$temp_4, na.rm = T),
                          temp_5 = mean(af_sig1$temp_5, na.rm = T),
                          temp_6 = mean(af_sig1$temp_6, na.rm = T),
                          temp_7 = mean(af_sig1$temp_7, na.rm = T),
                          
                          df_1 = mean(af_sig1$df_1, na.rm = T),
                          df_2 = mean(af_sig1$df_2, na.rm = T),
                          df_3 = mean(af_sig1$df_3, na.rm = T),
                          df_4 = mean(af_sig1$df_4, na.rm = T),
                          df_5 = mean(af_sig1$df_5, na.rm = T),
                          df_6 = mean(af_sig1$df_6, na.rm = T),
                          df_7 = mean(af_sig1$df_7, na.rm = T),
                          
                          idf_1 = mean(af_sig1$idf_1, na.rm = T),
                          idf_2 = mean(af_sig1$idf_2, na.rm = T),
                          idf_3 = mean(af_sig1$idf_3, na.rm = T),
                          idf_4 = mean(af_sig1$idf_4, na.rm = T),
                          idf_5 = mean(af_sig1$idf_5, na.rm = T),
                          idf_6 = mean(af_sig1$idf_6, na.rm = T),
                          idf_7 = mean(af_sig1$idf_7, na.rm = T),
                          
                          ied_1 = mean(af_sig1$ied_1, na.rm = T),
                          ied_2 = mean(af_sig1$ied_2, na.rm = T),
                          ied_3 = mean(af_sig1$ied_3, na.rm = T),
                          ied_4 = mean(af_sig1$ied_4, na.rm = T),
                          ied_5 = mean(af_sig1$ied_5, na.rm = T),
                          ied_6 = mean(af_sig1$ied_6, na.rm = T),
                          ied_7 = mean(af_sig1$ied_7, na.rm = T),
                          
                          week = median(af_sig1$week),
                          month = median(af_sig1$month),
                          station = 409900)
###
### Then, get predictions:
###

# The following code is sourced from Stack Overflow:
# https://stackoverflow.com/questions/30491545/predict-method-for-felm-from-lfe-package

predict.felm <- function(object, newdata, se.fit = FALSE,
                         interval = "none",
                         level = 0.95){
  if(missing(newdata)){
    stop("predict.felm requires newdata and predicts for all group effects = 0.")
  }
  
  tt <- terms(object)
  Terms <- delete.response(tt)
  attr(Terms, "intercept") <- 0
  
  m.mat <- model.matrix(Terms, data = newdata)
  m.coef <- as.numeric(object$coef)
  fit <- as.vector(m.mat %*% object$coef)
  fit <- data.frame(fit = fit)
  
  if(se.fit | interval != "none"){
    if(!is.null(object$clustervcv)){
      vcov_mat <- object$clustervcv
    } else if (!is.null(object$robustvcv)) {
      vcov_mat <- object$robustvcv
    } else if (!is.null(object$vcv)){
      vcov_mat <- object$vcv
    } else {
      stop("No vcv attached to felm object.")
    }
    se.fit_mat <- sqrt(diag(m.mat %*% vcov_mat %*% t(m.mat)))
  }
  if(interval == "confidence"){
    t_val <- qt((1 - level) / 2 + level, df = object$df.residual)
    fit$lwr <- fit$fit - t_val * se.fit_mat
    fit$upr <- fit$fit + t_val * se.fit_mat
  } else if (interval == "prediction"){
    stop("interval = \"prediction\" not yet implemented")
  }
  if(se.fit){
    return(list(fit=fit, se.fit=se.fit_mat))
  } else {
    return(fit)
  }
}

pred_IQ.w.u.l <- predict(IQ.w.u.l, new_data_iq, interval="confidence")
pred_IQ.w.u.q <- predict(IQ.w.u.q, new_data_iq, interval="confidence")
pred_IQ.w.u.c <- predict(IQ.w.u.c, new_data_iq, interval="confidence")

pred_IQ.m.u.l <- predict(IQ.m.u.l, new_data_iq, interval="confidence")
pred_IQ.m.u.q <- predict(IQ.m.u.q, new_data_iq, interval="confidence")
pred_IQ.m.u.c <- predict(IQ.m.u.c, new_data_iq, interval="confidence")

pred_IQ.w.c.l <- predict(IQ.w.c.l, new_data_iq, interval="confidence")
pred_IQ.w.c.q <- predict(IQ.w.c.q, new_data_iq, interval="confidence")
pred_IQ.w.c.c <- predict(IQ.w.c.c, new_data_iq, interval="confidence")

pred_IQ.m.c.l <- predict(IQ.m.c.l, new_data_iq, interval="confidence")
pred_IQ.m.c.q <- predict(IQ.m.c.q, new_data_iq, interval="confidence")
pred_IQ.m.c.c <- predict(IQ.m.c.c, new_data_iq, interval="confidence")

pred_IQ.w.i.l <- predict(IQ.w.i.l, new_data_iq, interval="confidence")
pred_IQ.w.i.q <- predict(IQ.w.i.q, new_data_iq, interval="confidence")
pred_IQ.w.i.c <- predict(IQ.w.i.c, new_data_iq, interval="confidence")

pred_IQ.m.i.l <- predict(IQ.m.i.l, new_data_iq, interval="confidence")
pred_IQ.m.i.q <- predict(IQ.m.i.q, new_data_iq, interval="confidence")
pred_IQ.m.i.c <- predict(IQ.m.i.c, new_data_iq, interval="confidence")

# Afghanistan:
pred_AF.w.u.l <- predict(AF.w.u.l, new_data_af, interval="confidence")
pred_AF.w.u.q <- predict(AF.w.u.q, new_data_af, interval="confidence")
pred_AF.w.u.c <- predict(AF.w.u.c, new_data_af, interval="confidence")

pred_AF.m.u.l <- predict(AF.m.u.l, new_data_af, interval="confidence")
pred_AF.m.u.q <- predict(AF.m.u.q, new_data_af, interval="confidence")
pred_AF.m.u.c <- predict(AF.m.u.c, new_data_af, interval="confidence")

pred_AF.w.i.l <- predict(AF.w.i.l, new_data_af, interval="confidence")
pred_AF.w.i.q <- predict(AF.w.i.q, new_data_af, interval="confidence")
pred_AF.w.i.c <- predict(AF.w.i.c, new_data_af, interval="confidence")

pred_AF.m.i.l <- predict(AF.m.i.l, new_data_af, interval="confidence")
pred_AF.m.i.q <- predict(AF.m.i.q, new_data_af, interval="confidence")
pred_AF.m.i.c <- predict(AF.m.i.c, new_data_af, interval="confidence")

###
### Then, build dataset from the predictions (for statistically significant models only:
###

preds_ols <- cbind(pred_IQ.w.u.q, 
                   pred_IQ.m.u.q, 
                   pred_AF.w.u.l, pred_AF.w.u.q, pred_AF.w.u.c, 
                   pred_AF.m.u.l, pred_AF.m.u.q, pred_AF.m.u.c)
colnames(preds_ols) <- c("model1_iq1", "model2_iq2", 
                         "model3_af1", "model4_af2", "model5_af3", 
                         "model6_af4", "model7_af5", "model8_af6")
# Create temperature range:
preds_ols$temperature <- temperature_range

# Take the means:
preds_ols$mean_iq <- rowMeans(preds_ols[, 1:2])
preds_ols$mean_af <- rowMeans(preds_ols[, 3:8])

# Then, distinguish predictions for temperatures outside of each country's observed values:
temperature_range_iq <- round(min(iq_sig1$temp)):round(max(iq_sig1$temp))
temperature_range_af <- round(min(af_sig1$temp)):round(max(af_sig1$temp))

preds_ols$model1_iq1_c <- 
  preds_ols$model2_iq2_c <- 
  preds_ols$model3_af1_c <- preds_ols$model4_af2_c <- preds_ols$model5_af3_c <- 
  preds_ols$model6_af4_c <- preds_ols$model7_af5_c <- preds_ols$model8_af6_c <- NA

preds_ols[which(preds_ols$temperature %in% temperature_range_iq), ]$model1_iq1_c <- preds_ols[which(preds_ols$temperature %in% temperature_range_iq), ]$model1_iq1
preds_ols[which(preds_ols$temperature %in% temperature_range_iq), ]$model2_iq2_c <- preds_ols[which(preds_ols$temperature %in% temperature_range_iq), ]$model2_iq2
preds_ols[which(preds_ols$temperature %in% temperature_range_af), ]$model3_af1_c <- preds_ols[which(preds_ols$temperature %in% temperature_range_af), ]$model3_af1
preds_ols[which(preds_ols$temperature %in% temperature_range_af), ]$model4_af2_c <- preds_ols[which(preds_ols$temperature %in% temperature_range_af), ]$model4_af2
preds_ols[which(preds_ols$temperature %in% temperature_range_af), ]$model5_af3_c <- preds_ols[which(preds_ols$temperature %in% temperature_range_af), ]$model5_af3
preds_ols[which(preds_ols$temperature %in% temperature_range_af), ]$model6_af4_c <- preds_ols[which(preds_ols$temperature %in% temperature_range_af), ]$model6_af4
preds_ols[which(preds_ols$temperature %in% temperature_range_af), ]$model7_af5_c <- preds_ols[which(preds_ols$temperature %in% temperature_range_af), ]$model7_af5
preds_ols[which(preds_ols$temperature %in% temperature_range_af), ]$model8_af6_c <- preds_ols[which(preds_ols$temperature %in% temperature_range_af), ]$model8_af6

preds_ols$mean_iq_c <- NA
preds_ols[which(preds_ols$temperature %in% temperature_range_iq), ]$mean_iq_c <- as.vector(preds_ols[which(preds_ols$temperature %in% temperature_range_iq), ]$mean_iq)
preds_ols$mean_af_c <- NA
preds_ols[which(preds_ols$temperature %in% temperature_range_af), ]$mean_af_c <- preds_ols[which(preds_ols$temperature %in% temperature_range_af), ]$mean_af

###
### Then, plot:
###

p_combined <- ggplot() + 
  geom_histogram(data = all_observed_temps[all_observed_temps$temperature %in% temperature_range, , drop=F], aes(x=temperature, y=(..ncount..)), bins = 100, size=.1, color="black", alpha=.4) +
  
  geom_line(data = preds_ols, aes(x = temperature, y = model1_iq1), color = "blue4", size = 0.15, linetype = "dashed") +
  geom_line(data = preds_ols, aes(x = temperature, y = model2_iq2), color = "blue4", size = 0.15, linetype = "dotted") +
  geom_line(data = preds_ols, aes(x = temperature, y = model3_af1), color = "deepskyblue2", size = 0.15, linetype = "dashed") +
  geom_line(data = preds_ols, aes(x = temperature, y = model4_af2), color = "deepskyblue2", size = 0.15, linetype = "dotted") +
  geom_line(data = preds_ols, aes(x = temperature, y = model5_af3), color = "deepskyblue2", size = 0.15, linetype = "dashed") +
  geom_line(data = preds_ols, aes(x = temperature, y = model6_af4), color = "deepskyblue2", size = 0.15, linetype = "dotted") +
  geom_line(data = preds_ols, aes(x = temperature, y = model7_af5), color = "deepskyblue2", size = 0.15, linetype = "dashed") +
  geom_line(data = preds_ols, aes(x = temperature, y = model8_af6), color = "deepskyblue2", size = 0.15, linetype = "dashed") +
  
  geom_line(data = preds_ols, aes(x = temperature, y = model1_iq1_c), color = "blue4", size = 0.5, linetype = "dashed") +
  geom_line(data = preds_ols, aes(x = temperature, y = model2_iq2_c), color = "blue4", size = 0.5, linetype = "dotted") +
  geom_line(data = preds_ols, aes(x = temperature, y = model3_af1_c), color = "deepskyblue2", size = 0.5, linetype = "dashed") +
  geom_line(data = preds_ols, aes(x = temperature, y = model4_af2_c), color = "deepskyblue2", size = 0.5, linetype = "dotted") +
  geom_line(data = preds_ols, aes(x = temperature, y = model5_af3_c), color = "deepskyblue2", size = 0.5, linetype = "dashed") +
  geom_line(data = preds_ols, aes(x = temperature, y = model6_af4_c), color = "deepskyblue2", size = 0.5, linetype = "dotted") +
  geom_line(data = preds_ols, aes(x = temperature, y = model7_af5_c), color = "deepskyblue2", size = 0.5, linetype = "dashed") +
  geom_line(data = preds_ols, aes(x = temperature, y = model8_af6_c), color = "deepskyblue2", size = 0.5, linetype = "dotted") +
  
  geom_line(data = preds_ols, aes(x = temperature, y = mean_iq_c), color = "blue4", size = 1) +
  geom_line(data = preds_ols, aes(x = temperature, y = mean_af_c), color = "deepskyblue2", size = 1) +
  
  geom_line(data = preds_ols, aes(x = temperature, y = mean_iq), color = "blue4", size = 0.5) +
  geom_line(data = preds_ols, aes(x = temperature, y = mean_af), color = "deepskyblue2", size = 0.5) +
  
  theme_bw() +
  geom_vline(xintercept = preds_ols[which(preds_ols$mean_iq == max(preds_ols$mean_iq)), ]$temperature, linetype=4, color = "blue4") +
  geom_vline(xintercept = preds_ols[which(preds_ols$mean_af == max(preds_ols$mean_af)), ]$temperature-0.25, linetype=4, color = "deepskyblue2") +
  
  ggtitle("") +
  ylab("Predicted Instances of (Least Constrained) Violent Attacks") + 
  xlab("Temperature (°F)") +
  xlim(-20,120) + 
  ylim(0,4) + 
  
  scale_x_continuous(breaks = seq(-20, 120, by = 5)) + 
  theme(
    axis.ticks = element_blank(), 
    text = element_text(size=18),
    axis.text.x = element_text(angle = 45))
p_combined

# Before generating predicted counts, check robustness to using an inverse hyperbolic sine
# transformation of the outcome variables instead (IHS: ln(x + \sqrt(x^2 + 1))):

iq_sig1$ihs_u <- log(iq_sig1$unconstrained + sqrt(iq_sig1$unconstrained^2 + 1))
iq_sig1$ihs_c <- log(iq_sig1$constrained + sqrt(iq_sig1$constrained^2 + 1))
iq_sig1$ihs_i <- log(iq_sig1$ied + sqrt(iq_sig1$ied^2 + 1))

af_sig1$ihs_u <- log(af_sig1$df + sqrt(af_sig1$df^2 + 1))
af_sig1$ihs_i <- log(af_sig1$ied + sqrt(af_sig1$ied^2 + 1))

# Iraq:
IQ.w.u.l_ihs <- felm(ihs_u ~ temp + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + ied + constrained +
                       temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 + 
                       constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                       unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                       idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                       ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                       week + governorate | 0 | 0, data = iq_sig1)
IQ.w.u.q_ihs <- felm(ihs_u ~ temp + temp2 + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + ied + constrained +
                       temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                       constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                       unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                       idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                       ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                       week + governorate | 0 | 0, data = iq_sig1)
IQ.w.u.c_ihs <- felm(ihs_u ~ temp + temp2 + temp3 + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + ied + constrained +
                       temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                       constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                       unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                       idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                       ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                       week + governorate | 0 | 0, data = iq_sig1)

IQ.w.c.l_ihs <- felm(ihs_c ~ temp + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + ied + unconstrained +
                       temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                       constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                       unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                       idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                       ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                       week + governorate | 0 | 0, data = iq_sig1)
IQ.w.c.q_ihs <- felm(ihs_c ~ temp + temp2 + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + ied + unconstrained +
                       temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                       constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                       unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                       idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                       ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                       week + governorate | 0 | 0, data = iq_sig1)
IQ.w.c.c_ihs <- felm(ihs_c ~ temp + temp2 + temp3 + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + ied + unconstrained +
                       temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                       constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                       unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                       idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                       ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                       week + governorate | 0 | 0, data = iq_sig1)

IQ.w.i.l_ihs <- felm(ihs_i ~ temp + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + constrained + unconstrained +
                       temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                       constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                       unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                       idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                       ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                       week + governorate | 0 | 0, data = iq_sig1)
IQ.w.i.q_ihs <- felm(ihs_i ~ temp + temp2 + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + constrained + unconstrained +
                       temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                       constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                       unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                       idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                       ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                       week + governorate | 0 | 0, data = iq_sig1)
IQ.w.i.c_ihs <- felm(ihs_i ~ temp + temp2 + temp3 + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + constrained + unconstrained +
                       temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                       constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                       unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                       idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                       ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                       week + governorate | 0 | 0, data = iq_sig1)

# Afghanistan:
AF.w.u.l_ihs <- felm(ihs_u ~ temp + prcp + wdsp + dewp + visib + mxspd + daylight + idf + ied +
                       temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                       df_1 + df_2 + df_3 + df_4 + df_5 + df_6 + df_7 +
                       idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                       ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                       week + station | 0 | month, data = af_sig1)
AF.w.u.q_ihs <- felm(ihs_u ~ temp + temp2 + prcp + wdsp + dewp + visib + mxspd + daylight + idf + ied +
                       temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                       df_1 + df_2 + df_3 + df_4 + df_5 + df_6 + df_7 +
                       idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                       ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                       week + station | 0 | month, data = af_sig1)
AF.w.u.c_ihs <- felm(ihs_u ~ temp + temp2 + temp3 + prcp + wdsp + dewp + visib + mxspd + daylight + idf + ied +
                       temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                       df_1 + df_2 + df_3 + df_4 + df_5 + df_6 + df_7 +
                       idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                       ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                       week + station | 0 | month, data = af_sig1)

AF.w.i.l_ihs <- felm(ihs_i ~ temp + prcp + wdsp + dewp + visib + mxspd + daylight + df + idf +
                       temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                       df_1 + df_2 + df_3 + df_4 + df_5 + df_6 + df_7 +
                       idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                       ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                       week + station | 0 | month, data = af_sig1)
AF.w.i.q_ihs <- felm(ihs_i ~ temp + temp2 + prcp + wdsp + dewp + visib + mxspd + daylight + df + idf +
                       temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                       df_1 + df_2 + df_3 + df_4 + df_5 + df_6 + df_7 +
                       idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                       ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                       week + station | 0 | month, data = af_sig1)
AF.w.i.c_ihs <- felm(ihs_i ~ temp + temp2 + temp3 + prcp + wdsp + dewp + visib + mxspd + daylight + df + idf +
                       temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                       df_1 + df_2 + df_3 + df_4 + df_5 + df_6 + df_7 +
                       idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                       ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                       week + station | 0 | month, data = af_sig1)

stargazer(IQ.w.u.l_ihs, IQ.w.u.q_ihs, IQ.w.u.c_ihs, 
          IQ.w.c.l_ihs, IQ.w.c.q_ihs, IQ.w.c.c_ihs, 
          IQ.w.i.l_ihs, IQ.w.i.q_ihs, IQ.w.i.c_ihs)

stargazer(AF.w.u.l_ihs, AF.w.u.q_ihs, AF.w.u.c_ihs, 
          AF.w.i.l_ihs, AF.w.i.q_ihs, AF.w.i.c_ihs)

###
### END.
###

# Next, generate predicted counts:

###
### IRAQ.
###

# Baghdad and Basrah:
# Unconstrained Violence:
# QP model results:

IQ.w.u.l.qp <- bayesglm(unconstrained ~ temp + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + ied + constrained +
                          temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                          constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                          unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                          idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                          ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                          as.factor(week) + as.factor(governorate), family = "quasipoisson", data = iq_sig1)
IQ.w.u.q.qp <- bayesglm(unconstrained ~ temp + temp2 + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + ied + constrained +
                          temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                          constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                          unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                          idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                          ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                          as.factor(week) + as.factor(governorate), family = "quasipoisson", data = iq_sig1)
IQ.w.u.c.qp <- bayesglm(unconstrained ~ temp + temp2 + temp3 + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + ied + constrained +
                          temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                          constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                          unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                          idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                          ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                          as.factor(week) + as.factor(governorate), family = "quasipoisson", data = iq_sig1)

IQ.m.u.l.qp <- bayesglm(unconstrained ~ temp + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + ied + constrained +
                          temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                          constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                          unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                          idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                          ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                          as.factor(month) + as.factor(governorate), family = "quasipoisson", data = iq_sig1)
IQ.m.u.q.qp <- bayesglm(unconstrained ~ temp + temp2 + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + ied + constrained +
                          temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                          constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                          unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                          idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                          ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                          as.factor(month) + as.factor(governorate), family = "quasipoisson", data = iq_sig1)
IQ.m.u.c.qp <- bayesglm(unconstrained ~ temp + temp2 + temp3 + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + ied + constrained +
                          temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                          constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                          unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                          idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                          ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                          as.factor(month) + as.factor(governorate), family = "quasipoisson", data = iq_sig1)

# Constrained Violence:
IQ.w.c.l.qp <- bayesglm(constrained ~ temp + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + ied + unconstrained +
                          temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                          constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                          unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                          idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                          ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                          as.factor(week) + as.factor(governorate), family = "quasipoisson", data = iq_sig1)
IQ.w.c.q.qp <- bayesglm(constrained ~ temp + temp2 + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + ied + unconstrained +
                          temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                          constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                          unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                          idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                          ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                          as.factor(week) + as.factor(governorate), family = "quasipoisson", data = iq_sig1)
IQ.w.c.c.qp <- bayesglm(constrained ~ temp + temp2 + temp3 + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + ied + unconstrained +
                          temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                          constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                          unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                          idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                          ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                          as.factor(week) + as.factor(governorate), family = "quasipoisson", data = iq_sig1)

IQ.m.c.l.qp <- bayesglm(constrained ~ temp + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + ied + unconstrained +
                          temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                          constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                          unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                          idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                          ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                          as.factor(month) + as.factor(governorate), family = "quasipoisson", data = iq_sig1)
IQ.m.c.q.qp <- bayesglm(constrained ~ temp + temp2 + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + ied + unconstrained +
                          temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                          constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                          unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                          idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                          ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                          as.factor(month) + as.factor(governorate), family = "quasipoisson", data = iq_sig1)
IQ.m.c.c.qp <- bayesglm(constrained ~ temp + temp2 + temp3 + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + ied + unconstrained +
                          temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                          constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                          unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                          idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                          ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                          as.factor(month) + as.factor(governorate), family = "quasipoisson", data = iq_sig1)

# IED:
IQ.w.i.l.qp <- bayesglm(ied ~ temp + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + constrained + unconstrained +
                          temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                          constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                          unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                          idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                          ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                          as.factor(week) + as.factor(governorate), family = "quasipoisson", data = iq_sig1)
IQ.w.i.q.qp <- bayesglm(ied ~ temp + temp2 + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + constrained + unconstrained +
                          temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                          constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                          unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                          idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                          ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                          as.factor(week) + as.factor(governorate), family = "quasipoisson", data = iq_sig1)
IQ.w.i.c.qp <- bayesglm(ied ~ temp + temp2 + temp3 + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + constrained + unconstrained +
                          temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                          constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                          unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                          idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                          ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                          as.factor(week) + as.factor(governorate), family = "quasipoisson", data = iq_sig1)

IQ.m.i.l.qp <- bayesglm(ied ~ temp + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + constrained + unconstrained +
                          temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                          constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                          unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                          idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                          ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                          as.factor(month) + as.factor(governorate), family = "quasipoisson", data = iq_sig1)
IQ.m.i.q.qp <- bayesglm(ied ~ temp + temp2 + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + constrained + unconstrained +
                          temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                          constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                          unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                          idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                          ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                          as.factor(month) + as.factor(governorate), family = "quasipoisson", data = iq_sig1)
IQ.m.i.c.qp <- bayesglm(ied ~ temp + temp2 + temp3 + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + constrained + unconstrained +
                          temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                          constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                          unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                          idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                          ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                          as.factor(month) + as.factor(governorate), family = "quasipoisson", data = iq_sig1)

###
### AFGHANISTAN.
###

AF.w.u.l.qp <- bayesglm(df ~ temp + prcp + wdsp + dewp + visib + mxspd + daylight + idf + ied +
                          temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                          df_1 + df_2 + df_3 + df_4 + df_5 + df_6 + df_7 +
                          idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                          ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                          as.factor(week) + as.factor(district), family = "quasipoisson", data = af_sig1)
AF.w.u.q.qp <- bayesglm(df ~ temp + temp2 + prcp + wdsp + dewp + visib + mxspd + daylight + idf + ied +
                          temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                          df_1 + df_2 + df_3 + df_4 + df_5 + df_6 + df_7 +
                          idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                          ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                          as.factor(week) + as.factor(district), family = "quasipoisson", data = af_sig1)
AF.w.u.c.qp <- bayesglm(df ~ temp + temp2 + temp3 + prcp + wdsp + dewp + visib + mxspd + daylight + idf + ied +
                          temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                          df_1 + df_2 + df_3 + df_4 + df_5 + df_6 + df_7 +
                          idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                          ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                          as.factor(week) + as.factor(district), family = "quasipoisson", data = af_sig1)

AF.m.u.l.qp <- bayesglm(df ~ temp + prcp + wdsp + dewp + visib + mxspd + daylight + idf + ied +
                          temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                          df_1 + df_2 + df_3 + df_4 + df_5 + df_6 + df_7 +
                          idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                          ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                          as.factor(month) + as.factor(district), family = "quasipoisson", data = af_sig1)
AF.m.u.q.qp <- bayesglm(df ~ temp + temp2 + prcp + wdsp + dewp + visib + mxspd + daylight + idf + ied +
                          temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                          df_1 + df_2 + df_3 + df_4 + df_5 + df_6 + df_7 +
                          idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                          ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                          as.factor(month) + as.factor(district), family = "quasipoisson", data = af_sig1)
AF.m.u.c.qp <- bayesglm(df ~ temp + temp2 + temp3 + prcp + wdsp + dewp + visib + mxspd + daylight + idf + ied +
                          temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                          df_1 + df_2 + df_3 + df_4 + df_5 + df_6 + df_7 +
                          idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                          ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                          as.factor(month) + as.factor(district), family = "quasipoisson", data = af_sig1)

AF.w.i.l.qp <- bayesglm(ied ~ temp + prcp + wdsp + dewp + visib + mxspd + daylight + df + idf +
                          temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                          df_1 + df_2 + df_3 + df_4 + df_5 + df_6 + df_7 +
                          idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                          ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                          as.factor(week) + as.factor(district), family = "quasipoisson", data = af_sig1)
AF.w.i.q.qp <- bayesglm(ied ~ temp + temp2 + prcp + wdsp + dewp + visib + mxspd + daylight + df + idf +
                          temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                          df_1 + df_2 + df_3 + df_4 + df_5 + df_6 + df_7 +
                          idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                          ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                          as.factor(week) + as.factor(district), family = "quasipoisson", data = af_sig1)
AF.w.i.c.qp <- bayesglm(ied ~ temp + temp2 + temp3 + prcp + wdsp + dewp + visib + mxspd + daylight + df + idf +
                          temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                          df_1 + df_2 + df_3 + df_4 + df_5 + df_6 + df_7 +
                          idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                          ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                          as.factor(week) + as.factor(district), family = "quasipoisson", data = af_sig1)

AF.m.i.l.qp <- bayesglm(ied ~ temp + prcp + wdsp + dewp + visib + mxspd + daylight + df + idf +
                          temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                          df_1 + df_2 + df_3 + df_4 + df_5 + df_6 + df_7 +
                          idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                          ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                          as.factor(month) + as.factor(district), family = "quasipoisson", data = af_sig1)
AF.m.i.q.qp <- bayesglm(ied ~ temp + temp2 + prcp + wdsp + dewp + visib + mxspd + daylight + df + idf +
                          temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                          df_1 + df_2 + df_3 + df_4 + df_5 + df_6 + df_7 +
                          idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                          ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                          as.factor(month) + as.factor(district), family = "quasipoisson", data = af_sig1)
AF.m.i.c.qp <- bayesglm(ied ~ temp + temp2 + temp3 + prcp + wdsp + dewp + visib + mxspd + daylight + df + idf +
                          temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                          df_1 + df_2 + df_3 + df_4 + df_5 + df_6 + df_7 +
                          idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                          ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                          as.factor(month) + as.factor(district), family = "quasipoisson", data = af_sig1)

# Next, generate predicted counts:

# First, construct a function for generating counts and creating output that can be plotted.
plot_function <- function(reg_output, X, dataset, type){
  set.seed(1)
  beta.sim <- mvrnorm(100, mu = coef(reg_output), Sigma = vcov(reg_output))
  paste <- model.matrix(reg_output)
  names1 <- colnames(paste)
  beta.sim <- as.data.frame(beta.sim)
  beta.sim <- beta.sim[, c(names1)]
  
  reps.df1 <- matrix(NA, length(unique(round(dataset$temp))), 100)
  for(i in 1:length(unique(round(dataset$temp)))){
    if(X %in% "linear"){
      paste[, 2] <- i + (min(round(dataset$temp))-1)
      x.beta <- paste %*% t(beta.sim)
      yhat <- exp(x.beta)
      reps.df1[i, ] <- apply(yhat, 2, mean, na.rm=T)
    } else if(X %in% "quadratic"){
      
      paste[, 2] <- i + (min(round(dataset$temp))-1)
      paste[, 3] <- (i + (min(round(dataset$temp))-1))^2
      x.beta <- paste %*% t(beta.sim)
      yhat <- exp(x.beta)

      reps.df1[i, ] <- apply(yhat, 2, mean, na.rm=T)
    } else if(X %in% "cubic"){
      paste[, 2] <- i + (min(round(dataset$temp))-1)
      paste[, 3] <- (i + (min(round(dataset$temp))-1))^2
      paste[, 4] <- (i + (min(round(dataset$temp))-1))^3
      x.beta <- paste %*% t(beta.sim)
      yhat <- exp(x.beta)
      reps.df1[i, ] <- apply(yhat, 2, mean, na.rm=T)
    }
  }
  
  plot_dataset <- as.data.frame(matrix(NA, length(unique(round(dataset$temp))),3))
  colnames(plot_dataset) <- c("means", "cil", "ciu")
  
  plot_dataset$means <- apply(reps.df1, 1 , mean)
  se <- apply(reps.df1, 1 , sd)
  plot_dataset$cil <- plot_dataset$means + qnorm(0.975)*se
  plot_dataset$ciu <- plot_dataset$means - qnorm(0.975)*se
  plot_dataset$temperature <- sort(unique(round(dataset$temp)))
  
  if(X %in% "linear"){
    plot_dataset$model <- "linear"
  } else if(X %in% "quadratic"){
    plot_dataset$model <- "quadratic"
  } else if(X %in% "cubic"){
    plot_dataset$model <- "cubic"
  }
  
  plot_dataset$max <- 0
  plot_dataset[which.max(plot_dataset$means), ]$max <- 1
  plot_dataset$attack_type <- type
  
  return(plot_dataset)
}

# Iraq:

p.IQ.w.u.q.qp <- plot_function(IQ.w.u.q.qp, "quadratic", iq_sig1, "unconstrained")
p.IQ.w.u.c.qp <- plot_function(IQ.w.u.c.qp, "cubic", iq_sig1, "unconstrained")

p.IQ.m.u.q.qp <- plot_function(IQ.m.u.q.qp, "quadratic", iq_sig1, "unconstrained")
p.IQ.m.u.c.qp <- plot_function(IQ.m.u.c.qp, "cubic", iq_sig1, "unconstrained")

p.IQ.w.c.l.qp <- plot_function(IQ.w.c.l.qp, "linear", iq_sig1, "constrained")
p.IQ.w.c.q.qp <- plot_function(IQ.w.c.q.qp, "quadratic", iq_sig1, "constrained")
p.IQ.w.c.c.qp <- plot_function(IQ.w.c.c.qp, "cubic", iq_sig1, "constrained")

p.IQ.m.c.l.qp <- plot_function(IQ.m.c.l.qp, "linear", iq_sig1, "constrained")
p.IQ.m.c.q.qp <- plot_function(IQ.m.c.q.qp, "quadratic", iq_sig1, "constrained")
p.IQ.m.c.c.qp <- plot_function(IQ.m.c.c.qp, "cubic", iq_sig1, "constrained")

p.IQ.w.i.q.qp <- plot_function(IQ.w.i.q.qp, "quadratic", iq_sig1, "ied")
p.IQ.w.i.c.qp <- plot_function(IQ.w.i.c.qp, "cubic", iq_sig1, "ied")

p.IQ.m.i.q.qp <- plot_function(IQ.m.i.q.qp, "quadratic", iq_sig1, "ied")
p.IQ.m.i.c.qp <- plot_function(IQ.m.i.c.qp, "cubic", iq_sig1, "ied")

# Afghanistan:

p.AF.w.u.l.qp <- plot_function(AF.w.u.l.qp, "linear", af_sig1, "unconstrained")
p.AF.w.u.q.qp <- plot_function(AF.w.u.q.qp, "quadratic", af_sig1, "unconstrained")
p.AF.w.u.c.qp <- plot_function(AF.w.u.c.qp, "cubic", af_sig1, "unconstrained")

p.AF.m.u.l.qp <- plot_function(AF.m.u.l.qp, "linear", af_sig1, "unconstrained")
p.AF.m.u.q.qp <- plot_function(AF.m.u.q.qp, "quadratic", af_sig1, "unconstrained")
p.AF.m.u.c.qp <- plot_function(AF.m.u.c.qp, "cubic", af_sig1, "unconstrained")

p.AF.w.i.l.qp <- plot_function(AF.w.i.l.qp, "linear", af_sig1, "ied")
p.AF.w.i.q.qp <- plot_function(AF.w.i.q.qp, "quadratic", af_sig1, "ied")
p.AF.w.i.c.qp <- plot_function(AF.w.i.c.qp, "cubic", af_sig1, "ied")

p.AF.m.i.l.qp <- plot_function(AF.m.i.l.qp, "linear", af_sig1, "ied")
p.AF.m.i.q.qp <- plot_function(AF.m.i.q.qp, "quadratic", af_sig1, "ied")
p.AF.m.i.c.qp <- plot_function(AF.m.i.c.qp, "cubic", af_sig1, "ied")

# Plot:

# Iraq:
# Least Constrained:
plot.IQ.m.u.l.qp <- ggplot() + 
  geom_line(data = p.IQ.w.u.q.qp, aes(y = means, x = temperature)) + 
  geom_ribbon(data = p.IQ.w.u.q.qp, aes(x = temperature, y = means, ymin = cil, ymax = ciu), alpha = .15, fill = "blue4") +
  geom_line(data = p.IQ.w.u.c.qp, aes(y = means, x = temperature)) + 
  geom_ribbon(data = p.IQ.w.u.c.qp, aes(x = temperature, y = means, ymin = cil, ymax = ciu), alpha = .15, fill = "blue4") +
  
  geom_line(data = p.IQ.m.u.q.qp, aes(y = means, x = temperature), color = "gray20", alpha = 0.5) + 
  geom_ribbon(data = p.IQ.m.u.q.qp, aes(x = temperature, y = means, ymin = cil, ymax = ciu), alpha = .15, fill = "gray20") +
  geom_line(data = p.IQ.m.u.c.qp, aes(y = means, x = temperature), color = "gray20", alpha = 0.5) +
  geom_ribbon(data = p.IQ.m.u.c.qp, aes(x = temperature, y = means, ymin = cil, ymax = ciu), alpha = .15, fill = "gray20") +
  
  theme_bw() +
  theme(text = element_text(size=26), legend.title=element_blank(),
        panel.background=element_blank(),
        panel.grid.minor=element_blank(), plot.background=element_blank()) +
  theme(legend.key = element_blank()) +
  scale_x_continuous(breaks = seq(min(round(af_sig1$temp)), max(round(af_sig1$temp)), by = 6)) + 
  theme(
    axis.ticks = element_blank(), 
    text = element_text(size=12),
    axis.text.x = element_text(angle = 45)) + 
  ylim(2, 10) +
  ylab("") + 
  xlab("") +
  ggtitle("Least Constrained")

# Most Constrained:
plot.IQ.m.c.l.qp <- ggplot() + 
  geom_line(data = p.IQ.w.c.l.qp, aes(y = means, x = temperature)) +
  geom_ribbon(data = p.IQ.w.c.l.qp, aes(x = temperature, y = means, ymin = cil, ymax = ciu), alpha = .15, fill = "green4") +
  geom_line(data = p.IQ.w.c.q.qp, aes(y = means, x = temperature)) + 
  geom_ribbon(data = p.IQ.w.c.q.qp, aes(x = temperature, y = means, ymin = cil, ymax = ciu), alpha = .15, fill = "green4") +
  geom_line(data = p.IQ.w.c.c.qp, aes(y = means, x = temperature)) + 
  geom_ribbon(data = p.IQ.w.c.c.qp, aes(x = temperature, y = means, ymin = cil, ymax = ciu), alpha = .15, fill = "green4") +
  
  geom_line(data = p.IQ.m.c.l.qp, aes(y = means, x = temperature), color = "gray20", alpha = 0.5) +
  geom_ribbon(data = p.IQ.m.c.l.qp, aes(x = temperature, y = means, ymin = cil, ymax = ciu), alpha = .15, fill = "gray20") +
  geom_line(data = p.IQ.m.c.q.qp, aes(y = means, x = temperature), color = "gray20", alpha = 0.5) + 
  geom_ribbon(data = p.IQ.m.c.q.qp, aes(x = temperature, y = means, ymin = cil, ymax = ciu), alpha = .15, fill = "gray20") +
  geom_line(data = p.IQ.m.c.c.qp, aes(y = means, x = temperature), color = "gray20", alpha = 0.5) +
  geom_ribbon(data = p.IQ.m.c.c.qp, aes(x = temperature, y = means, ymin = cil, ymax = ciu), alpha = .15, fill = "gray20") +
  
  theme_bw() +
  theme(text = element_text(size=26), legend.title=element_blank(),
        panel.background=element_blank(),
        panel.grid.minor=element_blank(), plot.background=element_blank()) +
  theme(legend.key = element_blank()) +
  scale_x_continuous(breaks = seq(min(round(af_sig1$temp)), max(round(af_sig1$temp)), by = 6)) + 
  theme(
    axis.ticks = element_blank(), 
    text = element_text(size=12),
    axis.text.x = element_text(angle = 45)) + 
  ylim(-0.5, 7.5) +
  ylab("") + 
  xlab("") +
  ggtitle("Most Constrained")

# IED:
plot.IQ.m.i.l.qp <- ggplot() + 
  geom_line(data = p.IQ.w.i.q.qp, aes(y = means, x = temperature)) + 
  geom_ribbon(data = p.IQ.w.i.q.qp, aes(x = temperature, y = means, ymin = cil, ymax = ciu), alpha = .15, fill = "firebrick") +
  geom_line(data = p.IQ.w.i.c.qp, aes(y = means, x = temperature)) + 
  geom_ribbon(data = p.IQ.w.i.c.qp, aes(x = temperature, y = means, ymin = cil, ymax = ciu), alpha = .15, fill = "firebrick") +
  
  geom_line(data = p.IQ.m.i.q.qp, aes(y = means, x = temperature), color = "gray20", alpha = 0.5) + 
  geom_ribbon(data = p.IQ.m.i.q.qp, aes(x = temperature, y = means, ymin = cil, ymax = ciu), alpha = .15, fill = "gray20") +
  geom_line(data = p.IQ.m.i.c.qp, aes(y = means, x = temperature), color = "gray20", alpha = 0.5) + 
  geom_ribbon(data = p.IQ.m.i.c.qp, aes(x = temperature, y = means, ymin = cil, ymax = ciu), alpha = .15, fill = "gray20") +
  
  theme_bw() +
  theme(text = element_text(size=26), legend.title=element_blank(),
        panel.background=element_blank(),
        panel.grid.minor=element_blank(), plot.background=element_blank()) +
  theme(legend.key = element_blank()) +
  scale_x_continuous(breaks = seq(min(round(af_sig1$temp)), max(round(af_sig1$temp)), by = 6)) + 
  theme(
    axis.ticks = element_blank(), 
    text = element_text(size=12),
    axis.text.x = element_text(angle = 45)) + 
  ylim(5, 13) +
  ylab("") + 
  xlab("") +
  ggtitle("IED Attacks")

# Afghanistan:
# DF:
plot.AF.m.u.l.qp <- ggplot() + 
  geom_line(data = p.AF.w.u.l.qp, aes(y = means, x = temperature)) +
  geom_ribbon(data = p.AF.w.u.l.qp, aes(x = temperature, y = means, ymin = cil, ymax = ciu), alpha = .15, fill = "blue4") +
  geom_line(data = p.AF.w.u.q.qp, aes(y = means, x = temperature)) + 
  geom_ribbon(data = p.AF.w.u.q.qp, aes(x = temperature, y = means, ymin = cil, ymax = ciu), alpha = .15, fill = "blue4") +
  geom_line(data = p.AF.w.u.c.qp, aes(y = means, x = temperature)) + 
  geom_ribbon(data = p.AF.w.u.c.qp, aes(x = temperature, y = means, ymin = cil, ymax = ciu), alpha = .15, fill = "blue4") +
  
  geom_line(data = p.AF.m.u.l.qp, aes(y = means, x = temperature), color = "gray20", alpha = 0.5) +
  geom_ribbon(data = p.AF.m.u.l.qp, aes(x = temperature, y = means, ymin = cil, ymax = ciu), alpha = .15, fill = "gray20") +
  geom_line(data = p.AF.m.u.q.qp, aes(y = means, x = temperature), color = "gray20", alpha = 0.5) + 
  geom_ribbon(data = p.AF.m.u.q.qp, aes(x = temperature, y = means, ymin = cil, ymax = ciu), alpha = .15, fill = "gray20") +
  geom_line(data = p.AF.m.u.c.qp, aes(y = means, x = temperature), color = "gray20", alpha = 0.5) + 
  geom_ribbon(data = p.AF.m.u.c.qp, aes(x = temperature, y = means, ymin = cil, ymax = ciu), alpha = .15, fill = "gray20") +
  
  theme_bw() +
  theme(text = element_text(size=26), legend.title=element_blank(),
        panel.background=element_blank(),
        panel.grid.minor=element_blank(), plot.background=element_blank()) +
  theme(legend.key = element_blank()) +
  scale_x_continuous(breaks = seq(min(round(af_sig1$temp)), max(round(af_sig1$temp)), by = 6)) + 
  theme(
    axis.ticks = element_blank(), 
    text = element_text(size=12),
    axis.text.x = element_text(angle = 45)) + 
  ylim(-0.5, 8) +
  ylab("") + 
  xlab("") +
  ggtitle("Direct Fire")

# IED:
plot.AF.m.i.l.qp <- ggplot() + 
  geom_line(data = p.AF.w.i.q.qp, aes(y = means, x = temperature)) + 
  geom_ribbon(data = p.AF.w.i.q.qp, aes(x = temperature, y = means, ymin = cil, ymax = ciu), alpha = .15, fill = "firebrick") +
  geom_line(data = p.AF.w.i.c.qp, aes(y = means, x = temperature)) + 
  geom_ribbon(data = p.AF.w.i.c.qp, aes(x = temperature, y = means, ymin = cil, ymax = ciu), alpha = .15, fill = "firebrick") +
  
  geom_line(data = p.AF.m.i.q.qp, aes(y = means, x = temperature), color = "gray20", alpha = 0.5) + 
  geom_ribbon(data = p.AF.m.i.q.qp, aes(x = temperature, y = means, ymin = cil, ymax = ciu), alpha = .15, fill = "gray20") +
  geom_line(data = p.AF.m.i.c.qp, aes(y = means, x = temperature), color = "gray20", alpha = 0.5) + 
  geom_ribbon(data = p.AF.m.i.c.qp, aes(x = temperature, y = means, ymin = cil, ymax = ciu), alpha = .15, fill = "gray20") +
  
  theme_bw() +
  theme(text = element_text(size=26), legend.title=element_blank(),
        panel.background=element_blank(),
        panel.grid.minor=element_blank(), plot.background=element_blank()) +
  theme(legend.key = element_blank()) +
  scale_x_continuous(breaks = seq(min(round(af_sig1$temp)), max(round(af_sig1$temp)), by = 6)) + 
  theme(
    axis.ticks = element_blank(), 
    text = element_text(size=12),
    axis.text.x = element_text(angle = 45)) + 
  ylim(-0.5, 8) +
  ylab("") + 
  xlab("") +
  ggtitle("IED Attacks")

# Combine Plots:

IQ.plots <- ggarrange(plot.IQ.m.u.l.qp + ylab("Iraq"), plot.IQ.m.c.l.qp, plot.IQ.m.i.l.qp, ncol = 3)
AF.plots <- ggarrange(plot.AF.m.u.l.qp + ylab("Afghanistan"), plot.AF.m.i.l.qp)
combined <- ggarrange(IQ.plots, AF.plots, nrow = 2)
combined <- annotate_figure(combined,
                            left = text_grob("Expected Attack Count", color = "black", rot = 90, size = 18),
                            bottom = text_grob("Temperature (°F)", color = "black", size = 18))

# Then, combine all figures:

# PRODUCES FIGURE 11 IN PAPER:
ggarrange(p_combined, combined)

###
### Baghdad Surveys
###

BG.m.l.b <- felm(vs_mnf ~ temp + age + work_wkhr + income.w + hhld_size + m_sig1_7 + prcp + wdsp + dewp + visib + mxspd + daylight +
                   p3_electr_binary + 
                   temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 |
                   education + block + date_month | 0 | block, data = bg_sur)
BG.m.q.b <- felm(vs_mnf ~ temp + temp2 + age + work_wkhr + income.w + hhld_size + m_sig1_7 + prcp + wdsp + dewp + visib + mxspd + daylight +
                   p3_electr_binary + 
                   temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 |
                   education + block + date_month | 0 | block, data = bg_sur)
BG.m.c.b <- felm(vs_mnf ~ temp + temp2 + temp3 + age + work_wkhr + income.w + hhld_size + m_sig1_7 + prcp + wdsp + dewp + visib + mxspd + daylight +
                   p3_electr_binary + 
                   temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 |
                   education + block + date_month | 0 | block, data = bg_sur)
BG.m.l.m <- felm(vs_mnf ~ temp + age + work_wkhr + income.w + hhld_size + m_sig1_7 + prcp + wdsp + dewp + visib + mxspd + daylight +
                   p3_electr_binary + 
                   temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 |
                   education + mahala + date_month | 0 | mahala, data = bg_sur)
BG.m.q.m <- felm(vs_mnf ~ temp + temp2 + age + work_wkhr + income.w + hhld_size + m_sig1_7 + prcp + wdsp + dewp + visib + mxspd + daylight +
                   p3_electr_binary + 
                   temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 |
                   education + mahala + date_month | 0 | mahala, data = bg_sur)
BG.m.c.m <- felm(vs_mnf ~ temp + temp2 + temp3 + age + work_wkhr + income.w + hhld_size + m_sig1_7 + prcp + wdsp + dewp + visib + mxspd + daylight +
                   p3_electr_binary + 
                   temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 |
                   education + mahala + date_month | 0 | mahala, data = bg_sur)
summary(BG.m.l.b)
summary(BG.m.q.b)
summary(BG.m.c.b)
summary(BG.m.l.m)
summary(BG.m.q.m)
summary(BG.m.c.m)

# How do you want to build this out?

stargazer(BG.m.q.b, BG.m.q.m)

BG.m.l.b.lg <- bayesglm(vs_mnf ~ temp + age + work_wkhr + income.w + hhld_size + m_sig1_7 + prcp + wdsp + dewp + visib + mxspd + as.factor(education) + as.factor(block) + as.factor(date_month), data = bg_sur, family = binomial(link = "logit"))
BG.m.q.b.lg <- bayesglm(vs_mnf ~ temp + temp2 + age + work_wkhr + income.w + hhld_size + m_sig1_7 + prcp + wdsp + dewp + visib + mxspd + as.factor(education) + as.factor(block) + as.factor(date_month), data = bg_sur, family = binomial(link = "logit"))
BG.m.c.b.lg <- bayesglm(vs_mnf ~ temp + temp2 + temp3 + age + work_wkhr + income.w + hhld_size + m_sig1_7 + prcp + wdsp + dewp + visib + mxspd + as.factor(education) + as.factor(block) + as.factor(date_month), data = bg_sur, family = binomial(link = "logit"))

BG.m.l.m.lg <- bayesglm(vs_mnf ~ temp + age + work_wkhr + income.w + hhld_size + m_sig1_7 + prcp + wdsp + dewp + visib + mxspd + as.factor(education) + as.factor(mahala) + as.factor(date_month), data = bg_sur, family = binomial(link = "logit"))
BG.m.q.m.lg <- bayesglm(vs_mnf ~ temp + temp2 + age + work_wkhr + income.w + hhld_size + m_sig1_7 + prcp + wdsp + dewp + visib + mxspd + as.factor(education) + as.factor(mahala) + as.factor(date_month), data = bg_sur, family = binomial(link = "logit"))
BG.m.c.m.lg <- bayesglm(vs_mnf ~ temp + temp2 + temp3 + age + work_wkhr + income.w + hhld_size + m_sig1_7 + prcp + wdsp + dewp + visib + mxspd + as.factor(education) + as.factor(mahala) + as.factor(date_month), data = bg_sur, family = binomial(link = "logit"))

#Next, repeat the analysis, using various income subsets:
bg_sur.s1.l <- bg_sur[bg_sur$income.w < 65000, ]
bg_sur.s2.l <- bg_sur[bg_sur$income.w < 70000, ]
bg_sur.s3.l <- bg_sur[bg_sur$income.w < 75000, ]
bg_sur.s4.l <- bg_sur[bg_sur$income.w < 80000, ]
bg_sur.s5.l <- bg_sur[bg_sur$income.w < 85000, ]
bg_sur.s1.u <- bg_sur[bg_sur$income.w > 65000, ]
bg_sur.s2.u <- bg_sur[bg_sur$income.w > 70000, ]
bg_sur.s3.u <- bg_sur[bg_sur$income.w > 75000, ]
bg_sur.s4.u <- bg_sur[bg_sur$income.w > 80000, ]
bg_sur.s5.u <- bg_sur[bg_sur$income.w > 85000, ]

BG.m.q.m.lg.s1.l <- bayesglm(vs_mnf ~ temp + temp2 + age + work_wkhr + income.w + hhld_size + m_sig1_7 + prcp + wdsp + dewp + visib + mxspd + as.factor(education) + as.factor(mahala) + as.factor(date_month), data = bg_sur.s1.l, family = binomial(link = "logit"))
BG.m.q.m.lg.s2.l <- bayesglm(vs_mnf ~ temp + temp2 + age + work_wkhr + income.w + hhld_size + m_sig1_7 + prcp + wdsp + dewp + visib + mxspd + as.factor(education) + as.factor(mahala) + as.factor(date_month), data = bg_sur.s2.l, family = binomial(link = "logit"))
BG.m.q.m.lg.s3.l <- bayesglm(vs_mnf ~ temp + temp2 + age + work_wkhr + income.w + hhld_size + m_sig1_7 + prcp + wdsp + dewp + visib + mxspd + as.factor(education) + as.factor(mahala) + as.factor(date_month), data = bg_sur.s3.l, family = binomial(link = "logit"))
BG.m.q.m.lg.s4.l <- bayesglm(vs_mnf ~ temp + temp2 + age + work_wkhr + income.w + hhld_size + m_sig1_7 + prcp + wdsp + dewp + visib + mxspd + as.factor(education) + as.factor(mahala) + as.factor(date_month), data = bg_sur.s4.l, family = binomial(link = "logit"))
BG.m.q.m.lg.s5.l <- bayesglm(vs_mnf ~ temp + temp2 + age + work_wkhr + income.w + hhld_size + m_sig1_7 + prcp + wdsp + dewp + visib + mxspd + as.factor(education) + as.factor(mahala) + as.factor(date_month), data = bg_sur.s5.l, family = binomial(link = "logit"))

BG.m.q.m.lg.s1.u <- bayesglm(vs_mnf ~ temp + temp2 + age + work_wkhr + income.w + hhld_size + m_sig1_7 + prcp + wdsp + dewp + visib + mxspd + as.factor(education) + as.factor(mahala) + as.factor(date_month), data = bg_sur.s1.u, family = binomial(link = "logit"))
BG.m.q.m.lg.s2.u <- bayesglm(vs_mnf ~ temp + temp2 + age + work_wkhr + income.w + hhld_size + m_sig1_7 + prcp + wdsp + dewp + visib + mxspd + as.factor(education) + as.factor(mahala) + as.factor(date_month), data = bg_sur.s2.u, family = binomial(link = "logit"))
BG.m.q.m.lg.s3.u <- bayesglm(vs_mnf ~ temp + temp2 + age + work_wkhr + income.w + hhld_size + m_sig1_7 + prcp + wdsp + dewp + visib + mxspd + as.factor(education) + as.factor(mahala) + as.factor(date_month), data = bg_sur.s3.u, family = binomial(link = "logit"))
BG.m.q.m.lg.s4.u <- bayesglm(vs_mnf ~ temp + temp2 + age + work_wkhr + income.w + hhld_size + m_sig1_7 + prcp + wdsp + dewp + visib + mxspd + as.factor(education) + as.factor(mahala) + as.factor(date_month), data = bg_sur.s4.u, family = binomial(link = "logit"))
BG.m.q.m.lg.s5.u <- bayesglm(vs_mnf ~ temp + temp2 + age + work_wkhr + income.w + hhld_size + m_sig1_7 + prcp + wdsp + dewp + visib + mxspd + as.factor(education) + as.factor(mahala) + as.factor(date_month), data = bg_sur.s5.u, family = binomial(link = "logit"))

plot_input_function <- function(x, y){
  set.seed(1)
  beta.sim <- mvrnorm(100, mu = coef(x), Sigma = vcov(x))
  paste <- model.matrix(x)
  names1 <- colnames(paste)
  beta.sim <- as.data.frame(beta.sim)
  beta.sim <- beta.sim[, c(names1)]
  
  reps <- matrix(NA, length(round(min(bg_sur$temp)):round(max(bg_sur$temp))), 100)
  for(i in 1:length(round(min(bg_sur$temp)):round(max(bg_sur$temp)))){
    if(y %in% "linear"){
      paste[, 2] <- i + (round(min(bg_sur$temp))-1)
      x.beta <- paste %*% t(beta.sim)
      yhat <- exp(x.beta) / (exp(x.beta) + 1)
      reps[i, ] <- apply(yhat, 2, mean, na.rm=T)
    } else if(y %in% "quadratic"){
      paste[, 2] <- i + (round(min(bg_sur$temp))-1)
      paste[, 3] <- i + (round(min(bg_sur$temp))-1)^2
      x.beta <- paste %*% t(beta.sim)
      yhat <- exp(x.beta) / (exp(x.beta) + 1)
      reps[i, ] <- apply(yhat, 2, mean, na.rm=T)
    } else if(y %in% "cubic"){
      paste[, 2] <- i + (round(min(bg_sur$temp))-1)
      paste[, 3] <- i + (round(min(bg_sur$temp))-1)^2
      paste[, 4] <- i + (round(min(bg_sur$temp))-1)^3
      x.beta <- paste %*% t(beta.sim)
      yhat <- exp(x.beta) / (exp(x.beta) + 1)
      reps[i, ] <- apply(yhat, 2, mean, na.rm=T)
    }
  }
  return(reps)
}

# First, plot temperature response results across functional forms:

# Next, given that the quadratic specification is the most significant, use that to 
# compare results across different income levels:

primary <- plot_input_function(BG.m.l.m.lg, "linear")
means.primary_l <- apply(primary, 1 , mean) 
se.primary <- apply(primary, 1 , sd)
cil.primary_l <- means.primary_l + qnorm(0.975)*se.primary
ciu.primary_l <- means.primary_l - qnorm(0.975)*se.primary
primary.dataset <- cbind(as.data.frame(means.primary_l), as.data.frame(cil.primary_l), as.data.frame(ciu.primary_l))

primary <- plot_input_function(BG.m.q.m.lg, "quadratic")
means.primary_q_m <- apply(primary, 1 , mean) 
se.primary <- apply(primary, 1 , sd)
cil.primary_q_m <- means.primary_q_m + qnorm(0.975)*se.primary
ciu.primary_q_m <- means.primary_q_m - qnorm(0.975)*se.primary
primary.dataset <- cbind(primary.dataset, as.data.frame(means.primary_q_m), as.data.frame(cil.primary_q_m), as.data.frame(ciu.primary_q_m))

primary <- plot_input_function(BG.m.q.b.lg, "quadratic")
means.primary_q_b <- apply(primary, 1 , mean) 
se.primary <- apply(primary, 1 , sd)
cil.primary_q_b <- means.primary_q_b + qnorm(0.975)*se.primary
ciu.primary_q_b <- means.primary_q_b - qnorm(0.975)*se.primary
primary.dataset <- cbind(primary.dataset, as.data.frame(means.primary_q_b), as.data.frame(cil.primary_q_b), as.data.frame(ciu.primary_q_b))

primary <- plot_input_function(BG.m.c.m.lg, "cubic")
means.primary_c <- apply(primary, 1 , mean) 
se.primary <- apply(primary, 1 , sd)
cil.primary_c <- means.primary_c + qnorm(0.975)*se.primary
ciu.primary_c <- means.primary_c - qnorm(0.975)*se.primary
primary.dataset <- cbind(primary.dataset, as.data.frame(means.primary_c), as.data.frame(cil.primary_c), as.data.frame(ciu.primary_c))
primary.dataset$temperature <- round(min(bg_sur$temp)):round(max(bg_sur$temp))

s1.l <- plot_input_function(BG.m.q.m.lg.s1.l, "quadratic")
means.s1.l <- apply(s1.l, 1 , mean) 
se.s1.l <- apply(s1.l, 1 , sd)
cil.s1.l <- means.s1.l + qnorm(0.975)*se.s1.l
ciu.s1.l <- means.s1.l - qnorm(0.975)*se.s1.l

s2.l <- plot_input_function(BG.m.q.m.lg.s2.l, "quadratic")
means.s2.l <- apply(s2.l, 1 , mean) 
se.s2.l <- apply(s2.l, 1 , sd)
cil.s2.l <- means.s2.l + qnorm(0.975)*se.s2.l
ciu.s2.l <- means.s2.l - qnorm(0.975)*se.s2.l

s3.l <- plot_input_function(BG.m.q.m.lg.s3.l, "quadratic")
means.s3.l <- apply(s3.l, 1 , mean) 
se.s3.l <- apply(s3.l, 1 , sd)
cil.s3.l <- means.s3.l + qnorm(0.975)*se.s3.l
ciu.s3.l <- means.s3.l - qnorm(0.975)*se.s3.l

s4.l <- plot_input_function(BG.m.q.m.lg.s4.l, "quadratic")
means.s4.l <- apply(s4.l, 1 , mean) 
se.s4.l <- apply(s4.l, 1 , sd)
cil.s4.l <- means.s4.l + qnorm(0.975)*se.s4.l
ciu.s4.l <- means.s4.l - qnorm(0.975)*se.s4.l

s5.l <- plot_input_function(BG.m.q.m.lg.s5.l, "quadratic")
means.s5.l <- apply(s5.l, 1 , mean) 
se.s5.l <- apply(s5.l, 1 , sd)
cil.s5.l <- means.s5.l + qnorm(0.975)*se.s5.l
ciu.s5.l <- means.s5.l - qnorm(0.975)*se.s5.l

s1.u <- plot_input_function(BG.m.q.m.lg.s1.u, "quadratic")
means.s1.u <- apply(s1.u, 1 , mean) 
se.s1.u <- apply(s1.u, 1 , sd)
cil.s1.u <- means.s1.u + qnorm(0.975)*se.s1.u
ciu.s1.u <- means.s1.u - qnorm(0.975)*se.s1.u

s2.u <- plot_input_function(BG.m.q.m.lg.s2.u, "quadratic")
means.s2.u <- apply(s2.u, 1 , mean) 
se.s2.u <- apply(s2.u, 1 , sd)
cil.s2.u <- means.s2.u + qnorm(0.975)*se.s2.u
ciu.s2.u <- means.s2.u - qnorm(0.975)*se.s2.u

s3.u <- plot_input_function(BG.m.q.m.lg.s3.u, "quadratic")
means.s3.u <- apply(s3.u, 1 , mean) 
se.s3.u <- apply(s3.u, 1 , sd)
cil.s3.u <- means.s3.u + qnorm(0.975)*se.s3.u
ciu.s3.u <- means.s3.u - qnorm(0.975)*se.s3.u

s4.u <- plot_input_function(BG.m.q.m.lg.s4.u, "quadratic")
means.s4.u <- apply(s4.u, 1 , mean) 
se.s4.u <- apply(s4.u, 1 , sd)
cil.s4.u <- means.s4.u + qnorm(0.975)*se.s4.u
ciu.s4.u <- means.s4.u - qnorm(0.975)*se.s4.u

s5.u <- plot_input_function(BG.m.q.m.lg.s5.u, "quadratic")
means.s5.u <- apply(s5.u, 1 , mean) 
se.s5.u <- apply(s5.u, 1 , sd)
cil.s5.u <- means.s5.u + qnorm(0.975)*se.s5.u
ciu.s5.u <- means.s5.u - qnorm(0.975)*se.s5.u

primary.dataset <- cbind(primary.dataset, 
                         as.data.frame(means.s1.l), as.data.frame(cil.s1.l), as.data.frame(ciu.s1.l),
                         as.data.frame(means.s2.l), as.data.frame(cil.s2.l), as.data.frame(ciu.s2.l),
                         as.data.frame(means.s3.l), as.data.frame(cil.s3.l), as.data.frame(ciu.s3.l),
                         as.data.frame(means.s4.l), as.data.frame(cil.s4.l), as.data.frame(ciu.s4.l),
                         as.data.frame(means.s5.l), as.data.frame(cil.s5.l), as.data.frame(ciu.s5.l),
                         as.data.frame(means.s1.u), as.data.frame(cil.s1.u), as.data.frame(ciu.s1.u), 
                         as.data.frame(means.s2.u), as.data.frame(cil.s2.u), as.data.frame(ciu.s2.u),
                         as.data.frame(means.s3.u), as.data.frame(cil.s3.u), as.data.frame(ciu.s3.u),
                         as.data.frame(means.s4.u), as.data.frame(cil.s4.u), as.data.frame(ciu.s4.u),
                         as.data.frame(means.s5.u), as.data.frame(cil.s5.u), as.data.frame(ciu.s5.u))

# Next, bound all negative predicted probabilities at 0 and 1:
primary.dataset <- primary.dataset %>% mutate(across(everything(), ~replace(., . < 0 , 0)))
primary.dataset <- primary.dataset %>% mutate(across(everything(), ~replace(., . > 1 , 1)))
# For plotting purposes, multiply all values by 100:
primary.dataset <- primary.dataset * 100

primary.dataset$temperature <-  round(min(bg_sur$temp)):round(max(bg_sur$temp))

distinct_survey_days <- distinct(bg_sur, date, .keep_all = TRUE)

p1 <- ggplot(data=primary.dataset) +

  geom_line(aes(x=temperature, y=means.primary_q_b), colour='blue4') +
  geom_ribbon(aes(x=temperature, ymin = cil.primary_q_b, ymax = ciu.primary_q_b), alpha = .2, fill = "blue4") +
  
  geom_line(aes(x=temperature, y=means.primary_q_m), colour='gray', alpha = 1) +
  geom_ribbon(aes(x=temperature, ymin = cil.primary_q_m, ymax = ciu.primary_q_m), alpha = .2, fill = "gray") +
  
theme_bw() +
  ggtitle("Full Iraqi Male Sample") +
  theme(text=element_text(size=14), 
        legend.title=element_blank(), 
        panel.background=element_blank(),
        panel.grid.minor=element_blank(), 
        plot.background=element_blank(), 
        axis.title.x=element_blank(),
        axis.title.y=element_blank()) +
  ylim(32, 100) +
  geom_vline(xintercept = which.is.max(means.primary_q_b)+(min(bg_sur$temp)-1), linetype=4) +
  xlim(40,105) + 
  scale_y_continuous(breaks = seq(32,100, by = 20), labels = scales::percent(seq(0.32,1, by = 0.20))) + 
  scale_x_continuous(breaks = seq(40,105, by = 5)) + 
  theme(
    axis.ticks = element_blank(), 
    text = element_text(size=14),
    axis.text.x = element_text(angle = 45))

p2 <- ggplot(data=primary.dataset) +
  geom_line(aes(x=temperature, y=means.s1.l, colour="< 65,000 IQD"), linetype="solid", alpha = 0.5) +
  geom_line(aes(x=temperature, y=means.s2.l, colour="< 70,000 IQD"), linetype="solid", alpha = 0.5) +
  geom_line(aes(x=temperature, y=means.s3.l, colour="< 75,000 IQD"), linetype="solid", alpha = 0.5) +
  geom_line(aes(x=temperature, y=means.s4.l, colour="< 80,000 IQD"), linetype="solid", alpha = 0.5) +
  geom_line(aes(x=temperature, y=means.s5.l, colour="< 85,000 IQD"), linetype="solid", alpha = 0.5) +
  
  geom_ribbon(aes(x=temperature, ymin = cil.s1.l, ymax = ciu.s1.l), alpha = .1, fill = "blue4") +
  geom_ribbon(aes(x=temperature, ymin = cil.s2.l, ymax = ciu.s2.l), alpha = .1, fill = "blue4") +
  geom_ribbon(aes(x=temperature, ymin = cil.s3.l, ymax = ciu.s3.l), alpha = .1, fill = "blue4") +
  geom_ribbon(aes(x=temperature, ymin = cil.s4.l, ymax = ciu.s4.l), alpha = .1, fill = "blue4") +
  geom_ribbon(aes(x=temperature, ymin = cil.s5.l, ymax = ciu.s5.l), alpha = .1, fill = "blue4") +
  
  ylim(32.5, 100) +
  
  theme_bw() +
  ggtitle("Highest Incomes Excluded") +
  scale_colour_manual("", 
                      breaks = c("< 65,000 IQD", "< 70,000 IQD", "< 75,000 IQD", "< 80,000 IQD", "< 85,000 IQD"),
                      values = c("firebrick", "green4", "blue4", "yellow4", "orange4")) +
  theme(text=element_text(size=14), 
        legend.title=element_blank(), 
        panel.background=element_blank(),
        panel.grid.minor=element_blank(), 
        plot.background=element_blank(), 
        axis.title.x=element_blank(),
        axis.title.y=element_blank())+
  geom_vline(xintercept = which.is.max(means.s3.l)+(min(bg_sur$temp)-1), linetype=4) +
  xlim(40,105) + 
  scale_y_continuous(breaks = seq(0,100, by = 20), labels = scales::percent(seq(0,1, by = 0.20))) + 
  scale_x_continuous(breaks = seq(40,105, by = 5)) + 
  theme(
    axis.ticks = element_blank(), 
    text = element_text(size=14),
    axis.text.x = element_text(angle = 45))

p3 <- ggplot(data=primary.dataset) +
  geom_histogram(data = distinct_survey_days, aes(x=temp, y=(..count..)), bins = 100, size=.1, color="black", alpha=.4) +
  
  geom_line(aes(x=temperature, y=means.s1.u, colour="> 65,000 IQD"), linetype="solid", alpha = 0.5) +
  geom_line(aes(x=temperature, y=means.s2.u, colour="> 70,000 IQD"), linetype="solid", alpha = 0.5) +
  geom_line(aes(x=temperature, y=means.s3.u, colour="> 75,000 IQD"), linetype="solid", alpha = 0.5) +
  geom_line(aes(x=temperature, y=means.s4.u, colour="> 80,000 IQD"), linetype="solid", alpha = 0.5) +
  geom_line(aes(x=temperature, y=means.s5.u, colour="> 85,000 IQD"), linetype="solid", alpha = 0.5) +
  
  geom_ribbon(aes(x=temperature, ymin = cil.s1.u, ymax = ciu.s1.u), alpha = .1, fill = "blue4") +
  geom_ribbon(aes(x=temperature, ymin = cil.s2.u, ymax = ciu.s2.u), alpha = .1, fill = "blue4") +
  geom_ribbon(aes(x=temperature, ymin = cil.s3.u, ymax = ciu.s3.u), alpha = .1, fill = "blue4") +
  geom_ribbon(aes(x=temperature, ymin = cil.s4.u, ymax = ciu.s4.u), alpha = .1, fill = "blue4") +
  geom_ribbon(aes(x=temperature, ymin = cil.s5.u, ymax = ciu.s5.u), alpha = .1, fill = "blue4") +
  ylim(0, 100) +
  theme_bw() +
  ggtitle("Highest Incomes Only") +
  scale_colour_manual("", 
                      breaks = c("> 65,000 IQD", "> 70,000 IQD", "> 75,000 IQD", "> 80,000 IQD", "> 85,000 IQD"),
                      values = c("firebrick", "green4", "blue4", "yellow4", "orange4")) +
  theme(text=element_text(size=14), 
        legend.title=element_blank(), 
        panel.background=element_blank(),
        panel.grid.minor=element_blank(), 
        plot.background=element_blank(), 
        axis.title.x=element_blank(),
        axis.title.y=element_blank())+
  geom_vline(xintercept = min(primary.dataset$temperature), linetype=4) +
  xlim(40,105) + 
  scale_y_continuous(breaks = seq(0,100, by = 20), labels = scales::percent(seq(0,1, by = 0.20))) + 
  scale_x_continuous(breaks = seq(40,105, by = 5)) + 
  theme(
    axis.ticks = element_blank(), 
    text = element_text(size=14),
    axis.text.x = element_text(angle = 45))

# PRODUCES FIGURE 14 IN PAPER:
plot.attitudes <- grid.arrange(p1, arrangeGrob(p2, p3, ncol=1), 
                               ncol=2, widths=c(1,1.2), bottom = textGrob("Temperature (°F)", gp=gpar(fontsize=18)), 
                               left = textGrob("Predicted Probability", rot = 90, gp=gpar(fontsize=18)), 
                               right = "", 
                               top="")

###
### END. 
###
